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1 Introduction 


2.1 Distribution on measurements 


In the biomagnetic inverse problem the main inter¬ 
est is the activation of a region of interest, i.e. the 
power dissipated in that region. The Bayesian power 
imaging method (BPI) provides a quantified probabil¬ 
ity that the activation of a region of interest is above 
a given threshold. This paper introduces the method 
and derives the equations used. The method is illus¬ 
trated in this paper using both experimental and sim¬ 
ulated data. Another paper [1] in this volume extends 
the method to compare task and control experiments. 


2 Methods 


Before we can begin to describe the method we need 
a few standard definitions, as used in [2]. Let m E 
R n be the vector of measurements at a single time 
instant, i.e. is the measurement of the ith channel, 
i = 1,... TV. Let Li(r) be the lead field of the ith 
detector. 

For simplicity we will use the same Hilbert space of 
currents as in [2], namely ^(Q) the square integrable 
field defined on the brain volume Q with the follow¬ 
ing inner product between currents 
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where uo(r) is a weighting distribution defined on the 
source space Q. Define the Gramm-Schmidt matrix 
P, using this inner product, by Pik = (pi , lj). 

This section continues by first deriving expressions 
for the a posteriori probability distribution on the 
measurements using the results from [2]. In turn this 
is used to derive an a posteriori probability distribu¬ 
tion on the activation of brain regions. It is this latter 
probability distribution that is mapped out in source 
space. 


The Bayesian method derived in [2] provides the a 
posteriori probability distribution of current distribu¬ 
tions j(r). The results of [2] are stated in terms of the 
projection of this a posteriori probability distribution 
onto one dimensional subspaces of Iq,(Q) by using 
a test current t(r). A special case (namely when the 
assumed prior current distribution j prior is identically 
zero) of the results derived in [2] are that the expected 
mean is u T (P + (D) _1 m and the expected variance 
is: 
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where u % = (£(r) , uo{r)Li{r)j . a 2 D is the covari¬ 
ance matrix of the assumed Gaussian noise and ( is a 
regularization parameter. In this Bayesian setting the 
regularization parameter ( is equal to a?//3 2 where 
f3 2 is the variance of the assumed prior probability 
distribution on currents j(r). 

In [2] the test current t was scanned across a set of 
voxels to produce an image. Here the approach is dif¬ 
ferent. By setting t(r) = oo(fr)Li(r) we can project 
the a posteriori distribution onto the expected mea¬ 
surement in the ith channel. So, the probability dis¬ 
tribution represents our best guess of the true value of 
the measurements given our prior knowledge of the 
measurement geometry and the data vector m. The 
expected mean and variance can be calculated using 
the equations above. The formulae can be simplified 
by noting that Uj = (uj(r)Lfr) , uo(r)Lj(r) ) — Pip 
i.e. the vector u is the ith column of P. It follows that 
the mean and variance for all channels can be calcu¬ 
lated in a single matrix equation. In this case the ex¬ 
pected mean is P(P + C,D)~ l m and the covariance 
matrix is given by 
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which simplifies to o?D(P + (D) l P. 



A few checks on the reasonableness of the above 
equations are possible by looking at extreme cases of 
the regularization parameter (. 

As ( —>► 0, the mean tends to m and the covariance 
tends to a 2 D. This is reasonable because, in this 
Bayesian setting, ( —>► 0 corresponds to having ei¬ 
ther perfect data (i.e. a 0) or no prior knowledge 
(i.e. (3 oo). In this case, the formulae show that 
the best estimate of the measurements is given by the 
data with covariance equal to the assumed covariance 
of the noise. 

As ( —>► oo, both the mean and the covariance matrix 
tend to zero. This is also reasonable since ( —>► oo 
corresponds to either worthless data (i.e. a oo) or 
perfect prior knowledge (i.e. (3 0). The formulae 

above indicate that the best estimate of the measure¬ 
ments is equal to our prior prejudices (in this case zero 
current density) with complete certainty. 

2.2 Distribution on activation estimates 

The above formulae can be used in the direct power 
estimation algorithm derived previously [3]. In this 
method an optimal estimation matrix Yq for the 0th 
voxel is derived and the expression rrfYom is used to 
estimate the activation of the 0th voxel. The expres¬ 
sion derived in [3] for Yq is 

(P 2 + D)Y e (P 2 +D) = PX e P (1) 

where Xq is a matrix which is related to the choice of 
region of interest. In this paper the region of interest is 
a single voxel in a source space so the expression for 
Xq is simply LqLq T where L k Q is a column of the gain 
matrix for the problem, i.e. is the measurement 
vector that would result from a dipole at position % 
oriented in direction e^. 

The matrix Y in Equation 1 is then given by; 
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where = (P 2 + D) 1 PLq. For this special case 
it is more efficient not to form the matrix Y but to 
calculate the activation directly using ; 
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AQ(m) = m T YQm = (3) 

k =o 

Here the method is put in a Bayesian context be¬ 
cause we now know the a posteriori distribution of 


measurement channels. Unfortunately the probabil¬ 
ity distribution of the quadratic form rPYQm is not 
known in a closed form. Analytic approximations 
have been derived to calculate the probability distri¬ 
bution of quadratic forms in Gaussian random vari¬ 
ables (e.g. [4]) but these methods are computationally 
expensive. In practice the quickest way of achieving a 
given level of accuracy for the probability is by Monte 
Carlo integration. First the measurements are normal¬ 
ized by factoring out the known covariance matrix: 

m * = -(P + CZ))- 1 /2pl/2 m ( 4 ) 
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Then the probability can be estimated by picking M 
Gaussian random vectors of length N with zero mean 
and unit covariance matrix, say z h , i = 1,... M, and 
computing the probability: 
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where Pq(x) is the one-sided probability of the esti¬ 
mate x and is the Heaviside step function. 

Maps of probability can then be computed by scan¬ 
ning 0 across a defined source space. 

3 Results 

The effectiveness of the algorithms has been tested in 
two simulated experiments and on real data. 

The first simulated experiment has a simple geometry 
(Figure 1). The head is modelled as a homogeneous 
conducting sphere of radius 8.9 cm with its centre at 
(0,0,—0.07cm). The source space is a 6cmx6cm 
square thin lamina consisting of 33x33 voxels in the 
plane £ = —0.01 cm with centre (0,0,-0.01 cm). The 
measurement instrument is a hexagonal array of 37 
second order axial gradiometers with baseline 5 cm 
with the lowest ’sensing’ coils in the plane £ = 4 cm. 

The sources for the first numerical study were 
two current dipoles with positions (0 cm, 0 cm) and 
(—2.25 cm,—2.25 cm) and oriented along the x and 
y axes respectively. Uncorrelated Gaussian noise is 
added to the pure signal so that the resulting signal to 
noise ratio was 9.7. 

In analyzing the experimental data there are many 
possible methods of determining noise levels etc. In 
this paper we have tried to keep the analysis choices 
simple in order to concentrate on the fundamentals 
of the method. With this in mind, all the examples 




Figure 1: A plan view of the experiment geometry. 
The grid depicts the source space voxels and the cir¬ 
cles represent the position of the measurement coils 
of the second order gradiometers. 

have the same choices for the parameters: the noise 
is assumed to be uncorrelated Gaussian noise of vari¬ 
ance a 2 , the regularization parameter is fixed to be 
0.1 x trace (P)/N and the number of trials, M, for 
the Monte-Carlo probability estimation method was 
1000. 

If the Monte-Carlo probability estimation method is 
applied directly to the raw power estimates, large ar¬ 
eas of the source space are found to be significant. 
This is because it is very unlikely for there to be zero 
power dissipated in each voxel. An easy calculation 
shows that the expected power for the 9th voxel is 
given by trace (Ye). The sensitivity to the (unknown) 
variance of the noise in the measurements can be min¬ 
imized by using the mean and standard deviation of 
the power estimates across the source space as a mea¬ 
sure of expected power. This is equivalent to the as¬ 
sumption that the desired signal produces significant 
power estimates on only a small proportion of the 
source space. In this paper the examples show prob¬ 
ability maps of the probability that the activation is 
more than three standard deviations above the mean 
activation. The resulting thresholded image is shown 
in Figure 2. The two sources are located accurately 
and with a very high significance. 

In the second example, a current dipole source was 
embedded in the innermost shell of a 3-shell model 
conducting sphere. The shell radii are 73 mm, 80 mm 
and 86 mm and the corresponding conductivities are 
0.33 Sm -1 , 0.0042 Sm -1 and 0.33 Sm -1 . A detector 
system consisting of 65 EEG electrodes was used to 
measure the scalp potential. The electrode positions 
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Figure 2: A probability map across the source space 
shown in Figure 1. The scale goes from white (repre¬ 
senting a probability of less than 0.1) to black (prob¬ 
ability greater than 0.9). 

and a representation of the potentials measured by the 
electrodes are shown in Figure 3. Also shown in Fig¬ 
ure 3 is the conducting shell used as the source space 
for the reconstruction. 



Figure 3: An EEG simulated dipolar dataset, together 
with the corresponding BPI map. The white blob that 
appears on the shell is the region where the activation 
is greater than the threshold (at the 99.9% signi cance 
level). 

The third and last example involves real MEG data 
from an evoked response study of face-processing [6]. 
The experimental instrument used is the Neuromag- 
122™ [5]. Human subjects were presented briefly 
with photographs of human faces. It is known that 
the early response to face images involves widespread 
activity in the posterior brain but there is limited ev- 



























idence for the precise distribution and timecourse of 
the neuronal sources. One suggestion is that there are 
three major areas of activity; in occipital cortex and 
both right and left ventral occipito-temporal cortex 
[7]. Strong occipital activity (starting about 100 ms 
after the stimulus) is expected to lead to concurrent 
activity in the two other regions with a stronger re¬ 
sponse in the right hemisphere [6]. 

The result shown in Figure 4 is derived from the re¬ 
sponse in an individual subject at a latency of 160 ms 
after the stimulus. At this time, there is a global 
peak in the total field power for this experiment. The 
source space was restricted to the surface of the cor¬ 
tex, as revealed by MRI and also shown in Figure 4. 
The findings of the BPI algorithm are consistent with 
earlier studies in that there is a significant response in 
both right and left cortices, with a stronger response 
in right hemisphere. Further work is needed to refine 
these predictions as the source space has a complex 
folded structure and we have not yet disentangled the 
contributions of superficial and deeper source regions 



Figure 4: The BPI map on the cortical source space 
used. The 122-channel helmet detector system is also 
shown. 


4 Discussion 

The main impetus behind the development of 
quadratic algorithms to estimate power rather than 
current density is that other modalities (e.g. fMRI, 
PET) naturally produce images of power dissipated. 
So if multi-modal comparisons are to be made then it 
is sensible to produce power images from MEG data. 
The BPI method has the additional advantage of be¬ 
ing able to produce spatial images of the variation of 
probability that are easy to interpret. 


In this paper the BPI method was shown to work ef¬ 
fectively for simulated MEG and EEG data. It was 
also demonstrated for experimental MEG data and 
produced results consistent with earlier analyses. 
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